nH0 <- 100
nH1 <- 200
phi_parametro <- 0.2
coeficientes <- function(phi) {
# βARMA: the model from Rocha and Cribari-Neto (2009, 2017)
# is obtained by setting coefs$d = 0 and d = FALSE and error.scale = 1 (predictive scale)
list(alpha = 0, beta = NULL, phi = phi, theta = NULL, d = 0, nu = 20)
}
barma.sim <- function(n, phi, seed, y.start = NULL){
BARFIMA.sim(
n = n,
coefs = coeficientes(phi),
y.start = y.start,
error.scale = 1,
complete = F,
seed = seed
)
}
barma.phi_estimado <- function(yt, alpha = 0, nu = 20, phi = 0.1){
BARFIMA.fit(
yt = yt,
start = list(alpha = alpha, nu = nu, phi = phi),
p = 1,
d = FALSE,
error.scale = 1,
report = F
)$coefficients["phi"][[1]]
}
barma.residuos <- function(yt, phi_estimado){
BARFIMA.extract(
yt = yt,
coefs = coeficientes(phi_estimado),
llk = F,
info = F,
error.scale = 1
)$error
}# Função para forçar a execução de `BARFIMA.extract` até que ela funcione
so_quero_que_funcione <- function(func, debug = T, max = 1000, ...) {
FINALMENTE <- F
contador <- 0
while (!FINALMENTE) {
contador <- contador + 1
if (debug) print(paste("Tentativa:", contador))
resultado <- tryCatch({
func(...)
}, error = function(err) {
if (contador >= max) {
stop("Deu ruim, número máximo de tentativas atingido")
}
if (any(grepl("\\.btsr\\.extract\\(", deparse(err$trace$call)))) {
# Ignora erro de extração de resíduos
return(NULL)
}
stop(err)
})
if (!is.null(resultado)) {
FINALMENTE <- T
}
}
return(resultado)
}ewma_qcc <- function(amostra_inicial, lambda, amostra_subsequente) {
ew <- ewma(amostra_inicial, lambda = lambda, nsigmas = 3, plot = F, newdata = amostra_subsequente)
registros <- (nH0 + 1):(nH0 + nH1)
list(
ewma = as.numeric(ew$y[registros]),
LI = ew$limits[, 1][registros],
LS = ew$limits[, 2][registros]
)
}gera_amostras_subsequentes <- function() {
lambda <- 0.2
amostra_controle <- barma.sim(
n = nH0,
phi = phi_parametro,
seed = 1337
)
ultima_observacao <- amostra_controle[nH0]
phi_estimado <- barma.phi_estimado(amostra_controle)
residuos_controle <- barma.residuos(amostra_controle, phi_estimado)
data.frame(
novo_phi = seq(0, 0.8, by = 0.1)
) %>%
rowwise() %>%
mutate(
observacao = list(1:nH1),
dados = list(barma.sim(
n = nH1,
phi = novo_phi,
seed = 1337,
y.start = ultima_observacao
)),
residuos = list(
barma.residuos(dados, phi_estimado)
),
qcc = list(ewma_qcc(residuos_controle, lambda, residuos)),
ewma = list(qcc$ewma),
LI = list(qcc$LI),
LS = list(qcc$LS),
fora_de_controle = list(ewma < LI | ewma > LS),
total_fora_de_controle = sum(fora_de_controle),
fracao_fora_de_controle = total_fora_de_controle / nH1
)
}
amostras_subsequentes <- so_quero_que_funcione(gera_amostras_subsequentes)## [1] "Tentativa: 1"amostras_subsequentes %>%
group_by(novo_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
.groups = "drop"
)ggplotly(
amostras_subsequentes %>%
unnest(
cols = c(novo_phi, observacao, ewma, LI, LS)
) %>%
ggplot() +
geom_line(aes(x = observacao, y = LI, color = "Limite Inferior"), linetype = "dashed") +
geom_line(aes(x = observacao, y = LS, color = "Limite Superior"), linetype = "dashed") +
geom_line(aes(x = observacao, y = ewma, color = "EWMA")) +
geom_point(data = . %>% filter(ewma < LI | ewma > LS),
aes(x = observacao, y = ewma, color = "Fora de controle"), size = 1.5, shape = 4) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
"Limite Superior" = adjustcolor("red", alpha.f = 0.5),
"EWMA" = adjustcolor("black", alpha.f = 0.8),
"Fora de controle" = adjustcolor("red", alpha.f = 0.4)
)
) +
labs(
title = "EWMA para resíduos βAR(1), com λ = 0.2",
) +
theme.base +
theme(
legend.position = "bottom",
legend.title = element_blank()
) +
facet_wrap(~novo_phi)
) %>%
plotly.baseamostras_monte_carlo_gerador <- function(){
numero_de_execucoes <- 200
expand.grid(
k = 1:numero_de_execucoes,
novo_phi = seq(0.2, 0.6, by = 0.1),
lambda = seq(0.1, 0.4, by = 0.1)
) %>%
mutate(
id = row_number()
) %>%
rowwise() %>%
mutate(
amostra_controle = list(barma.sim(
n = nH0,
phi = phi_parametro,
seed = id
)),
ultima_observacao = amostra_controle[nH0],
phi_estimado = barma.phi_estimado(amostra_controle),
residuos_controle = list(
barma.residuos(amostra_controle, phi_estimado)
),
dados = list(barma.sim(
n = nH1,
phi = novo_phi,
y.start = ultima_observacao,
seed = id + 1337E4
)),
residuos = list(
barma.residuos(dados, phi_estimado)
),
qcc = list(ewma_qcc(residuos_controle, lambda, residuos)),
ewma = list(qcc$ewma),
LI = list(qcc$LI),
LS = list(qcc$LS),
fora_de_controle = list(ewma < LI | ewma > LS),
total_fora_de_controle = sum(fora_de_controle, na.rm = TRUE),
fracao_fora_de_controle = total_fora_de_controle / nH1
)
}
if (file.exists("cache_simulacao_ewma.RData")) {
amostras_monte_carlo <- readRDS("cache_simulacao_ewma.RData")
} else {
amostras_monte_carlo <- so_quero_que_funcione(amostras_monte_carlo_gerador)
saveRDS(amostras_monte_carlo, file = "cache_simulacao_ewma.RData")
}amostras_monte_carlo_resumo <- amostras_monte_carlo %>%
group_by(lambda, novo_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
amostras_monte_carlo_resumoggplotly(
amostras_monte_carlo_resumo %>%
ggplot(aes(x = novo_phi, y = mean, color = factor(lambda))) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
color = "Valores de λ",
title = "Fração de pontos fora de controle por Φ e λ"
) +
theme.base
) %>%
plotly.baseggplotly(
amostras_monte_carlo %>%
ggplot(aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(lambda))) +
geom_boxplot() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
fill = "Valores de λ",
title = "Fração de pontos fora de controle por Φ e λ"
) +
theme.base
) %>%
plotly.base %>%
layout(boxmode = 'group')